Discrete Breathers in One-Dimensional Diatomic Granular Crystals 
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We report the experimental observation of discrete breathers in a one-dimensional diatomic gran- 
ular crystal composed of compressed elastic beads that interact via Hertzian contact. We first 
characterize their effective linear spectrum both theoretically and experimentally. We then illus- 
trate theoretically and numerically the modulational instability of the lower edge of the optical band. 
This leads to the dynamical formation of long-lived breather structures, whose families of solutions 
we compute throughout the linear spectral gap. Finally, we observe experimentally such localized 
breathing modes with quantitative characteristics that agree with our numerical results. 
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Introduction. Intrinsic localized modes (ILMs), or dis- 
crete breathers (DBs), have been a central theme in nu- 
merous theoretical and experimental investigations dur- 
ing the past two decades [l|, Q • Their original theoret- 
ical proposal in settings such as anharmonic nonlinear 
lattices [3j and the rigorous proof of their existence un- 
der fairly general conditions [4| motivated studies of such 
modes in a diverse host of applications, including charge- 
transfer solids [5], antiferromagnets [6[, superconducting 
Josephson junctions 0, photonic crystals (jl, biopoly- 
mers [9(, micromechanical cantilever arrays [lOj, Bosc- 
Einstein condensates [llj ]. and more. 

Granular crystals, consisting of closely packed ensem- 
bles of elastically interacting particles have also recently 
drawn considerable attention. This broad interest has 
arisen from their nonlinear contact dynamics and the 
tunability of their dynamic response to encompass linear, 
weakly nonlinear, and strongly nonlinear regimes 12, 13j |. 
Such flexibility makes them ideal not only as toy models 
for probing the physics of granular materials but also for 
the implementation of many engineering applications, in- 
cluding shock and energy absorbing layers [lj] , actuating 
devices [l5j], and sound scramblers Only recently 

have nonlinear localized modes begun to be explored 
in granular crystals. Previous studies have focused on 
metastable breathers in acoustic vacuum [lJljthe obser- 
vation of localized oscillations near a defect HI, , and 
one-dimensional (ID) diatomic crystals restricted to lin- 
ear dynamics due to welded sphere contacts [2(| • Under- 
standing and controlling localization in granular systems 
might lead to new energy harvesting/conversion devices 
and acoustic filters. 

In this Letter, we investigate the existence, stability, 
and dynamics of DBs in a compressed ID diatomic granu- 
lar crystal using experiments, theory, and numerical sim- 
ulations. We first detail our experimental setup and the- 
oretical model. We then analyze the system's dynamics 
in the linear regime, show how a modulational instabil- 
ity generates DBs in the weakly nonlinear regime, and 



finally provide experimental evidence of their existence. 

Experimental setup. We assemble a ID diatomic gran- 
ular crystal by alternating aluminum spheres (Acraball, 
6061-T6 type, with radius R a = 9.525 mm, mass m a — 
9.75 g, elastic modulus E a = 73.5 GPa, and Poisson ratio 
v a = 0.33) and stainless steel spheres (McMaster-Carr, 
316 type, R b = R a , m b = 28.84 g, E b = 193 GPa, and 
v b = 0.3). The reported values of E aib and u a<b are stan- 
dard specifications 2l| ; we discuss the precise character- 
ization of the effective elastic properties of our system 
below. We hold the spheres in place using four polycar- 
bonate restraining bars and guide plates. At one end 
of the crystal, we apply a precompressive force using a 
lever-mass system. We position an 0-1 tool steel plate, 
which we mount on a 1018 steel angle bracket at the 
other end of the crystal as a "wall" . We drive dynamical 
perturbations using a piezoelectric actuator, which we 
fit on the steel plate. We visualize the evolution of the 
force-time history of the propagating excitations using 
periodically-placed calibrated piezo sensors that we em- 
bed inside selected particles (preserving the inertia and 
the bulk stiffness of the original bead (la , [lij ] ) . We mea- 
sure the static load using a calibrated strain gauge cell 
that we place in contact with the lever arm and the last 
bead of the crystal. 

Theoretical model. We model a ID diatomic crystal of 
N spheres as a chain of nonlinear oscillators (l^ |: 
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where [Y]+ denotes the positive part of Y, Ui is the 
displacement of the ith sphere (where i S {1, • • • , N}) 
around the static equilibrium, the masses are m dd = fn a 
and m eve n = fnbi and the coefficient A depend on the 
exponent p and the geometry/material properties of ad- 
jacent beads. The exponent p = 3/2 yields the Hertz 
potential law between adjacent spheres (22j. In this case, 

A = (f ^ + V-l^y\-k + ^y^, and one ob- 
tains a static overlap of So = (Fq/A) 2 ^ 3 under a static 
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load Fq 13, 13| ■ We compute the linear dispersion curve 
of our system from the linearization of Eq. [1] For di- 
atomic crystals, this curve contains two branches (acous- 
tic and optical) 0. At the edge of the first Brillouin 
zone — i.e., at k = where a = R a + Rb — Sq is the 
equilibrium distance between two adjacent beads — the 
linear spectrum possesses a gap between the upper cutoff 
u>x = y^2K2 /M of the acoustic branch and the lower cut- 
off u>2 = \j2Kilm of the optical one; the linear stiffness 
is K 2 = |A 2 / 3 F^ 3 , and we define M — max {m a ,raf,} 
and m = min {m a , raj,}. The upper cutoff frequency of 
the optical band is located at W3 = y/2K2(l/m + 1/M). 
In Table HI we summarize K2, A, and the three cutoff 
frequencies, which we estimate using standard specifica- 
tions [2l| and compute using a static load of Fq = 20 N. 

Linear spectrum. We experimentally characterize the 
linear (phonon) spectrum of a diatomic crystal [24J (N — 
81 and Fq = 20 N) by applying low-amplitude (peak at 
approximately 10 mN), broadband (2 — 18 kHz frequency 
width) and uniform noise for 800 ms. We measure the 
dynamical forces using a sensor located inside the 14th 
particle, plus that from the driving voltage and the ac- 
tuator sensitivity. We then compute the power spectral 
density (PSD) [25| of the force-sensor, normalize it to 
the PSD of the driving force, and average the ratio over 
8 acquisitions to obtain the transfer function shown in 
Fig. Q] This spectrum clearly shows forbidden bands (i.e., 
gaps) and allowed bands bounded by cutoff frequencies. 
These frequencies match half of the transfer function's 
low-frequency level, which is determined as the average 
level in the 2 — 4 kHz range. We summarize these fre- 
quencies in Table [U Matching these frequencies to the 
theoretical formulas above provides an opportunity to 
probe the beads' effective parameters K2 and A shown in 
Table U (errorbars indicate the standard deviations from 
the three frequencies measurements) . We find that all of 
the cutoff frequencies show a systematic upshift of about 
9% compared to the predictions from standard specifica- 
tions. We identify four possible explanations for such a 
systematic bias: (i) the uncertainty in the standard val- 
ues of material parameters [2l| ; (ii) non-Hookean elastic 
dynamics might lead to slight shift on the nonlinear expo- 
nent p and accordingly a large deviation in the coefficient 
A [22|; (iii) imperfect surface smoothness might also in- 
duce fluctuations in p and hence A [26| ; and (iv) dissipa- 
tive mechanisms, such as viscoelasticity and solid friction, 
can induce stiffening of the interaction potential between 
particles 
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We also test a shorter (N = 15) crys- 
tal, which showed a higher low-frequency level and lower 
linear stiffness K2 — 14.80 N//im, in agreement with less 
dissipation in a shorter crystal. 

Modulational Instability and DBs. We now consider 
the weakly nonlinear dynamics of the granular crystal. If 
the displacements have small amplitudes relative to those 
due to precompression, we do a power series expansion 
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FIG. 1: [Color online] Top panel: Experimental setup. Bot- 
tom panel: Experimental phonon spectrum of the 81-bead 
steel-aluminum diatomic crystal. The horizontal line is half 
of the low frequency average level and vertical lines indicate 
the /£ xp cutoff frequencies given in Table HI 





fi [kHz] 
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h [kHz] 


K 2 [N//H 


A [N//im 3/2 ] 


th. 


4.71 


8.10 


9.37 


12.63 


5.46 


exp. 


5.11 


8.83 


10.22 


14.95 + 0.10 


7.04 + 0.07 


diff. 


+8.5% 


+9.0% 


+9.1% 


+18.4% 


+28.8% 



TABLE I: Predicted (from standard specifications [21|) ver- 
sus measured cutoff frequencies, linear stiffness K2, and coef- 
ficient A under a static precompression of Fq = 20 N. 



of the forces (up to quartic displacement terms) to yield 
the K2 — K3 — K4 model: 
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where K 3 = -|A 4 / 3 F 1/3 and K 4 
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Equa- 
tion {3J constitutes a diatomic variant of the Fermi- 
Pasta-Ulam (FPU) nonlinear oscillator chain [28]. Be- 

> |, the nonlinearity induces an asymmet- 
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ric localized mode in the gap of the linear spectrum 
(i.e., one obtains a gap soliton). This arises from the 
modulational instability (MI) of the optical lower cutoff 
phonon mode [29| . In order to verify this prediction, we 
solve Eq. ([T]) numerically using A exp (see Table HJ and 
the optical lower cutoff mode as initial condition. This 
mode corresponds to the crystal vibration in which the 
light masses oscillate with frequency f^ 9 and the heavy 
masses are at rest. In order to trigger the instability 
of this mode, we choose an oscillation amplitude of the 
light masses that corresponds to an 11.25 N dynamical 
peak force. As shown in Fig. [DJa), we observe the MI 
and the resulting generation of a localized mode with 
frequency /& ~ 7.95 kHz at t ~ 8 ms. In order to ob- 
serve the generation of DBs under conditions relevant to 
our experimental setup, wc also run simulations in which 
the displacement of the first (actuator-driven) aluminum 
bead is described by a 30 ms long square windowed sine 
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FIG. 2: [Color online] (a),(b): Spatiotemporal evolution of 
the forces [N] . (a) Manifestation of the MI of the optical lower 
cutoff mode, (b) Generation of a DB under conditions rele- 
vant to our experimental setup, (c), (d) PSDs of the forces 
at particles i = 36 and i = 32, respectively, for the DB simu- 
lation shown in panel (b). Dashed lines indicate the driving 
frequency / ac t = f^ P , and the arrows indicate the frequency 
component associated with the generated DB. 



FIG. 3: [Color online] Bifurcation diagram of the parameter 
continuation of the DB solutions, (a) Maximal force of the 
wave versus frequency ft, (along with profile insets at two val- 
ues of /;,). (b) Maximal deviation of Floquet multipliers from 
the unit circle, indicating the instability growth strength. The 
right inset shows a typical multiplier picture, and the left one 
shows the connection between the strong (real multiplier) in- 
stability and the change in sign of dE/dfb- 



with variable oscillation amplitude B and frequency at 
the optical lower cutoff / act = f^ 9 - In Fig. [2fb), we 
show an example of the spatiotemporal evolution of the 
forces when B = 0.061<5o- In this example, the maxi- 
mum dynamic force acting on the beads over the first 
10 cycles of the excitation is about 6.5 N ~ 0.325i*o- 
We thus anticipate a weakly nonlinear response that is 
well described by the K2 — K3 — K4 theory. Initially, a 
small amplitude excitation is generated and transmitted. 
During its transmission, the light masses oscillate out of 
phase, and the heavy ones are practically at rest. At 
t ~ 22ms, the excitation experiences MI, which yields a 
DB which, for these initial conditions, is localized near 
bead 37. This nonlinear solution exists even after the 
actuator is turned off at t — 30 ms. The PSDs of the 
force at particles i = 36 [see Fig. HJc)] and i = 32 [see 
Fig. [2jd)] reveal the presence of a frequency component 
in the gap at fb ~ 8.14 kHz < f^ v . Moreover, it is clear 
that the dominant frequency is fb at the center of the 
localized mode and that / act dominates 4 particles away. 

Exact solutions and stability of DBs. We apply New- 
ton's method (see [H and references therein) with free 
boundary conditions to numerically obtain, with high 
precision, the above dynamically generated DB wave- 
forms as exact (time-periodic) solutions. We then study 
their linear stability and frequency dependence (within 
the spectral gap). Continuing this solution within the 
gap [i.e., for / e (/T xp ,/2 Xp )] starting from the opti- 
cal cutoff mode allows us to trace the entire family of 
DB solutions. In Fig. [3fa), we show the maximum force 
max(i 7 i), which is the experimentally observable param- 
eter of the DB solution, as a function of the DB fre- 
quency fb- As fb — > f^ P , max(Fj) — > F and the DBs 
broaden and finally merge with the linear optical lower 
cutoff mode. In the insets of Fig. [3£a), we show exam- 
ples of these solutions with frequencies fbi = 8.35 kHz 



and fb2 = 8.75 kHz. To examine the stability of the DB 
solutions, we compute their Floquet multipliers Xj [2]. 
If \Xj\ = 1 for all j, then the DB is linearly stable. In 
Fig. [21b), we show the stability diagram for the family 
of DB solutions and the corresponding locations of Flo- 
quet multipliers in the complex plane for the DB with 
fbi = 8.83 kHz. Strictly speaking, the DB is stable only 
for fb ~ /2 XP - Otherwise, the DB family exhibits os- 
cillatory instabilities @, [l9|. However, the deviations of 
the unstable eigenvalues from the unit circle are bounded 
above by 0.08, and numerical integration of the DBs up 
to times 100T (where T is their period) reveals their ro- 
bustness. Importantly, we also find that DB solutions 
exhibit a strong instability due to a pair of real multipli- 
ers when f b S (8.45 kHz, 8.67 kHz). As seen in Fig.^b), 
this instability is connected with the turning points of the 
energy of the DB as a function of its frequency (these oc- 
cur when dE/dfb = 0). Similar features have also been 
observed in diatomic Klein-Gordon chains [30j. 

Experimental observation of DBs. Here, we excite the 
81-bead diatomic crystal with a higher-amplitude signal 
(relative to the linear-spectrum experiments). The actu- 
ator is driven by a 30 ms long square windowed sine volt- 
age. We probe a range of driving frequencies (near the 
lower optical cutoff frequency - see f% xp in Table HJ and 
amplitudes around the values expected to create DBs. 
We place force sensors in particles 2, 4, 7, 12, and 14. In 
Fig. 21 we show experimental evidence of a DB. For the 
case shown in Fig. 2J the driving frequency is 8.94 kHz 
and the peak force that we measure near the actuator 
is 12.20 N ~ 0.61F . Figure |4jb) shows the force versus 
time at particle 14, and the corresponding PSDs, shown 
in Figs. (Hd) and UJe) , demonstrate the existence of a 
second mode (whose frequency differs from the driving 
frequency), which we indicate with an arrow in the fig- 
ure. This mode occurs at /^ xp — 8.35 kHz, which is in the 
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FIG. 4: [Color online] Experimental observation of a DB at 
jcxp ^ g kjj z ( a ) Force at particle i = 2, (b) force at 
particle i = 14, (c) PSD at i = 2 for t > 1 ms, (d) PSD at 
i — 14 for 1 < i < 29 ms (while the actuator is on), and (e) 
PSD at i = 14 for t > 31 ms (after the actuator is switched 
off). Vertical lines in (a,b) mark the times when the actuator 
is switched on and off. Vertical lines in (c,d,e) indicate the 
driving frequency (8.94 kHz > f^ p ; see Table [I}, and Sf is 
the frequency resolution. Arrows in (d,e) indicate a DB that 
lies inside the band gap (/i Xp < /™ p < /| xp ; see Table IJ). 



band gap and yields a robust DB according to the previ- 
ous linear stability analysis. As with the PSD at particle 
2, which we show in Fig.[4jc), the PSDs at particles 4 and 
7 do not exhibit additional modes in the gap. The PSD 
of the signal from sensor 12 reveals the presence of this 
mode, but the amplitude is smaller than that from the 
sensor at i = 14, indicating that the center of the DB is 
located further inside the bulk of the crystal. Addition- 
ally, as predicted by simulations, the DB appears to be 
long-lived, as shown for instance from the different decay 
rates in Figs. HJa) and0Jb) after we switch off actuator. 
In Fig. H](e) , we estimate the PSD of the tail, illustrating 
that the DB maintains its prominence while the mode 
at the actuator frequency has experienced a decrease in 
PSD amplitude by two orders of magnitude. 

Conclusions. We have characterized the dynamics of 
compressed ID diatomic granular crystals using exper- 
iments, numerical simulations, and theoretical analysis. 
We found satisfactory agreement between experiments 
and theory for the linearized spectrum of steel-aluminum 
crystals. We also explored theoretically the formation 
of DBs via MI, which, in turn, led us to systematically 
trace them numerically and to observe them experimen- 
tally. Our results provide a first step towards achieving 
a deeper understanding and classifying intrinsic modes 
in ID granular crystals, and pave the way for their man- 
ifestation and dynamics in 2D and 3D settings, leading 
ultimately towards their potential exploitation in energy- 
harvesting applications. 
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